Gravitational waves from the Papaloizou-Pringle instability in black hole-torus 

systems 
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Black hole (BH)-torus systems are promising candidates for the central engine of gamma-ray 
bursts (GRBs), and also possible outcomes of the collapse of supermassive stars to supermassive 
black holes (SMBHs). By three-dimensional general relativistic numerical simulations, we show that 
an m = 1 nonaxisymmetric instability grows for a wide range of self-gravitating tori orbiting BHs. 
The resulting nonaxisymmetric structure persists for a timescale much longer than the dynamical 
one, becoming a strong emitter of large amplitude, quasiperiodic gravitational waves. Our results 
indicate that both, the central engine of GRBs and newly formed SMBHs, can be strong gravitational 
wave sources observable by forthcoming ground-based and spacecraft detectors. 

PACS numbers: 04.25.D-, 04.30.-w, 04.40.Dg 



Introduction. — Black hole (BH)-torus systems are 
common in the universe. The central region of active 
galactic nuclei is believed to consist of a supermassive 
black hole (SMBH) of mass A/bh ~ IO^-IO^OMo sur- 
rounded by a torus. Such systems may form through the 
collapse of supermassive stars (SMS) [H, 0. Mergers of 
neutron star (NS) binaries and BH-NS binaries often re- 
sult in a BH and a torus Q. Such systems can also be 
produced at the end of the life of massive stars [4] . It has 
been suggested that the merger and collapsar scenarios 
are linked to short- and long-duration gamma-ray bursts 
(GRBs), respectively [HI, Thus, BH-torus systems 
may become high-energy astrophysical objects. How- 
ever, their formation and evolution have not yet been ob- 
served as such sites are opaque to electromagnetic waves 
(EWs) due to their intrinsic high density and temper- 
ature. Gravitational waves (GWs), however, are much 
more transparent than EWs with respect to absorption 
and scattering with matter. If BH-torus systems emit- 
ted detectable GWs, it would be possible to explore their 
formation and evolution, along with the prevailing hy- 
potheses that associate them to GRB engines. This may 
become possible in the near future thanks to the array 
of ground-based (LIGO, VIRGO, LCGT) and spacecraft 
(LISA) detectors of GWs 0. 

BH-torus systems are known to be subject to the ax- 
isymmetric runaway instability and the nonaxisym- 
metric Papaloizou-Pringle instability (PPI) . Recently, 
general relativistic (GR) simulations have shown that the 
runaway instability does not have a significant impact 
on the dynamics even if the torus self-gravity is taken 
into account [loj . However, nonaxisymmetric features 
ma y s till play a crucial role as shown by linear analy- 
sis |11 | . Exploring the nonlinear growth and saturation of 
the PPI in BH-torus systems requires three-dimensional 
(3D) simulations [l2|. Moreover, since GR plays an im- 
portant role in the case of self-gravitating tori [13], 3D 



GR numerical simulations are mandatory. In particu- 
lar, long-term simulations are essential to follow the BH- 
torus systems after the PPI saturation. By contrast, the 
recent 3D GR simulations of BH-torus by Korobkin et 
al. only followed the linear growth phase of the PPI [l3| . 

In this Letter we report results of a sample of such 
simulations which show that the PPI sets in for a wide 
range of BH-torus systems, and that the resulting nonax- 
isymmetric structure is maintained after the saturation 
of the instability. As a result, such systems can be strong 
sources of detectable GWs for the detectors that will be in 
operation in the near future. We use units of c = G = 1 , 
and the perturbation is expressed by the expansion of 
exp{— i(cLii — rrnp)}. 

Method and initial models. — The simulations are 
performed using a fixed mesh refinement (FMR) and 
Message Passing Interface version of the SACRA code 



(see [lj| for details regarding the implementation of 



the FMR algorithm and the formulation and numerical 
schemes for solving Einstein's and hydrodynamics equa- 
tions). The FMR structure consists of 6 refinement levels 
for which the grid resolution changes as Axi-i — 2Axi 
{I — 2-6), with a finest level of Axi=e — 0.05AfBH- 
Each domain is composed of {2N + l,2N + l,N + 1) 
Cartesian grid zones (x, y, z) which cover the interval 
[—NAxi : N Axi] for the x— and y— directions and 
[0 : A^Aa;;] along z, with A^ = 98. The outer boundary is 
located at 156.8A/bh which is in the local wave zone (see 
Table H]). The BH horizon is covered by ~ 10 grid points 
in the finest refinement level. For assessing the validity 
of the numerical results, simulations with different values 
of Axi^g; and size of each domain were performed, 
which confirmed that convergence is achieved with our 
chosen values. We use a standard F-law equation of state 
(EOS) P — iV — l)pe, where P is the pressure, p the 
rest-mass density, e the specific internal energy, and F 
the adiabatic index. 
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TABLE I: List of key quantities for the initial data: specific angular momentum profile (value of 
j'/Mbh at the inner edge), torus-to-BH mass ratio, TZ = Mtor/A'/BH, position of the inner edge, 
center (location of the maximum density), and outer edge on the equatorial plane, rin, Vc, and rout, 
orbital period at r = rc, torb, and simulation time in units of torb- The last two columns show the 
growth rate lm{ujm) of the PPI and the Fourier amplitude of the density perturbation for the m — 1 
mode at the saturation, 5i = \Di\/Do. Here, 57c is the angular velocity at r = rc. 
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Compact equilibria for a BH-torus system are pre- 
pared in the puncture framework and used as initial con- 
ditions [ill- We adopt the T = 4/3 polytropic EOS, 



P — up , to mimic the EOS of a degenerate relativistic 
electron gas or radiation dominated gas. The polytropic 
constant n is chosen such that the torus-to-BH mass ra- 
tios are in the range 7^ = Mtor/MeH = 0.06-0.1. We 
specify the angular momentum profiles, j — where 
j = {l + e + P/ p)u^ and = u'^/u*, denoting the fluid 
four velocity. Table H] lists key quantities of the four dif- 
ferent BH-torus systems considered. The BH is assumed 
to be non-rotating. We investigate both j ^constant 
(C models) , and non-constant profiles (NC models) with 
j oc ri~^/^, to assess their infiuence. The selected mass 
ratios are partially motivated by the result of simulations 
of gravitational collapse to a BH and torus 0, [l^ . The 
numerical model of [l6| may be regarded as a candidate 
of the central engine of GRBs, and that of as a model 
for the formation of a SMBH through the collapse of a 
rotating SMS. 

Our numerical simulations are performed for a time 
long enough, i.e. t ^ 20-40 torb, to include the growth, 
saturation, and post-saturation phases of the PPI, for 
which the growth timescale is of order torb (see Table I). 

Results. — Figure [T] shows the evolution of the maxi- 
mum baryon rest-mass density of the torus, the baryon 
rest-mass fiux measured on the apparent horizon (AH), 
Mr^i-AHJ the total baryon rest-mass located outside 
the AH for each model. Figure ^ shows snapshots of the 
rest-mass density distribution on the equatorial plane af- 
ter the PPI saturation for models CI and NCI. Figure [3] 
shows the time evolution of the density perturbation am- 
plitudes for the TO = 1-3 modes of model CI, a typical 
power-law index of the angular velocity, and the outgo- 

(2 2) 

ing component of the complex Weyl scalar ' with 
I = m = 2 mode for models CI and NCI, where tret is 
the retarded time. 

These figures reveal that the PPI is present in all mod- 
els except NC06. The instability grows exponentially 
with time, with the m = 1 mode being the fastest grow- 
ing mode (cf. Fig. |3Ua)). Figure Ufa) indicates that dur- 
ing the linear growth phase, the maximum density is ap- 
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FIG. 1: Evolution of (a) the maximum rest-mass density, (b) 
the mass flux on the BH horizon, and (c) the rest-mass located 
outside the apparent horizon (AH) for each model. The max- 
imum density is evaluated outside the AH and rest-mass is 
normalized by its initial value. 



proximately constant, but after several orbits, nonlinear 
growth sets in and Pmax increases. After saturation, pmax 
decreases and eventually approaches a new quasiequilib- 
rium value. For a given angular momentum profile, the 
PPI growth timescale is shorter the larger the mass ratio 
TZ. For a given mass ratio, the growth timescale for NC 
models is longer than for C models. 

As shown in Fig. [IJb), the mass fiux increases expo- 
nentially with the PPI growth. It reaches its peak value 
at mode growth saturation and then settles to an ap- 
proximately non-zero constant value. These features are 
qualitatively identical for the PPI unstable models. The 
mass flux for NC06 is negligible throughout the evolution. 
These features are reflected in Fig.[IJc): Around the sat- 
uration time of the PPI nonlinear growth, the torus mass 
steeply decreases, and then the torus relaxes to a station- 
ary accretion state. For NC models, the mass flux at the 
saturation time is always lower than for C models, which 
indicates the dependence of the accretion history on the 
profile of j . In Fig. [5] we illustrate the nonaxisymmetric 
morphological features present in the tori of models CI 
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FIG. 2: Snapshots of the rest-mass density distribution on 
the equatorial plane for models CI at t = 15.11 torh (left) 
and NCI at t = 20.02 ta-h (right). Vectors show the velocity 
fields, and black filled circles around the centers show the 
region inside the AHs. 
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FIG. 3: Evolution of (a) the m = 1-3 mode amplitudes for 
model CI, (b) the power-law index of the angular velocity 
profile for all models, and (c) the outgoing part of the complex 
Weyl scalar Re(*f '^') for models CI (green) and NCI (blue). 



and NCI once the stationary accretion phase is reached. 
The presence of the m — 1 mode is apparent. These re- 
sults imply that more massive tori and/or & j — const 
profile favor the appearance of the PPI with respect to 
less massive tori and/or a j ^ const profile, for which 
the torus may be PP-stable (e.g. model NC06). In addi- 
tion, no sign of the runaway instability is found for any 
models. 

Figure[3l^a) plots the evolution of the Fourier modes for 
the density perturbation of model CI, for which we define 
Dm = / pe^^'^fcPx and Sm = \Dm\/Do. This shows 
that the m = 1 mode is the fastest growing mode and 
that saturation occurs at t ^ 10 torh- We determine the 
PPI growth rate through a fit of the numerical data, and 
find it in a wide range ~ 5-25% of the angular velocity 
at the saturation time (see Table H]) . As expected from 
Fig. [TJ massive and/or j — const models exhibit larger 
growth rates (in agreement with ) and no evidence of 
the instability is found for model NC06. 

The PPI activates the transport of angular momentum 
outwards through the corotation point resonance 11| . 
and therefore, redistributes angular momentum. Fig- 
ure Efb) shows the evolution of a typical power-law in- 
dex of the angular velocity profile, where we fit f2(r) = 
l/27r / n{r,TT/2,(p) dip by oc r'' for r G [^'im'^out] on the 
equatorial plane. In Newtonian gravity q — —2 for j = 
const and q =^ —1.5 for Keplerian motion. It is expected 
that q would approach —1.5 (Keplerian limit) from a 
lower initial value < —1.5. Note that in GR j = const law 
does not correspond to q = —2 and also q ^ —1.5 for the 
Keplerian limit. Figure [31[b) shows that after the onset 
of the PPI (for all models except for NC06), q increases 
and then saturates, refiecting the outward transport of 
angular momentum. The increase of q indicates that tori 
with a positive value of djjdr are stable against the PPI, 
as expected from linear analysis llj and from previous 
nonlinear simulations 13. All the features found here 



suggest that the PPI plays an important role for the BH- 
torus system dynamics unless the torus mass is small and 
the angular momentum profile is close to the Keplerian 
profile. 

A point worth stressing is that the m = 1 nonax- 
isymmetric structure survives with an appreciable am- 
plitude after the saturation of the PPI nonlinear growth, 
as shown in Fig. [5] The likely reason is that as the com- 
mon center of mass for the BH-torus system should re- 
main at the origin and the BH center is not located at 
the center of mass, then the torus should be deformed. 
The persistent m = 1 structure leads to the emission of 
quasiperiodic GWs with a large amplitude. Note that the 
dominant mode of GWs is / = m = 2 because the situa- 
tion is similar to a test particle orbiting around a BH. An 
outgoing component of the Weyl scalar extracted in the 
local wave zone, plotted in Fig. [31Jc) for model CI (green 
curve), shows that the waveforms refiect the evolution 
of the system: A burst at irot/iorb ~ 10 is emitted by 
the PPI nonlinear growth and saturation, and subsequent 
quasiperiodic waves for irct/^orb <^ 10 are emitted by the 
long-lived m = 1 structure. This feature is universal for 
all the C models. While the GW for the NCI model 
(Fig. Eljc), blue curve) does not show a prominent burst, 
they do show the quasiperiodic emission induced by the 
long-lived m = 1 structure. The striking modulation in 
the GW signal found only in j ^ const models may be due 
to the variability in the maximum rest-mass density and 
the g-index for the angular velocity after mode growth 
saturation (between lOiorb <t < BOtorb for model NCI). 
Therefore, it may be associated with the excitation of 
p-modes (pressure-inertial modes) in the torus. These p- 
mode oscillations are gradually damped, as well as the 
modulation, for t > SOtorb- The emission of quasiperi- 
odic GWs is expected to continue as far as the disk has 
a significant amount of mass. Such time duration can be 
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FIG. 4: GW spectra for all the models for Mbh = 10M@ and 
D = 100 Mpc, and Mbh = 10*^ Af© and D = 10 Gpc. The 
cross symbols show the peak amplitudes found in the numer- 
ical simulations and the circles show hypothetical amplitudes 
inferred from accretion timescales (see main text for details). 



estimated from the mass accretion rate, and the resid- 
ual disk mass (Fig. (rfb) and (c)). We find the accretion 
timescales are iacc ~ -^^r>rAH/-^r=rAH ~ x 10^ Mbh- 
In Fig. |4] we show with crosses the peak effective am- 
plitude hcs of GWs together with the design sensitivities 
of ground-based and spacecraft GW detectors [7] . Here 



^eff 



i{\h, 



\h^\^}/2)^/^D/MBH with 



and D 



being the Fourier component of two polarization modes 
and the distance from the source. Since our focus is on 
the GW emission after the PPI growth saturation, we dis- 
card the waveform before the saturation for the C models, 
e.g. for irct/^oi-b ^ 10 for CI, when calculating hcS- From 
the accretion timescale tacc, we may expect that the GW 
emission will continue for A'cydo ^ ^acc/^GW ~ 100 where 
tew is the typical oscillation period of GWs. Hence, the 
peak amplitudes may be enhanced by A'cydo 10, val- 
ues plotted with circles in Fig. SI Because we use units 
of Mbh = 1 in our simulations, the numerical results can 
be rescaled for arbitrary BH mass. We choose the hy- 
pothetical values for the collapsar model of GRB central 
engine as (A/bh, D) = (IOMq, 100 Mpc) and for the hy- 
pothetical formation scenario of a SMBH through SMS 
collapse as (Mbh,£') = (IO^Mq, 10 Gpc) in Fig.H 

For the stellar-mass BH, h^g could be ~ 10~^^ at the 
peak with the frequency 100-200 Hz, which corresponds 
to the Kepler frequency of rorb ~ 10A/bh- The peak 
amplitudes agree with the order estimation, 

2M^=lV^ r- 

/^pcak y-^ ''cycle 



D 

^ V in-21 \ / 7^ \ / Mbh \ /lOOMpc 



N, 



cycle 

100 



where M„j=i is the mass of matter which contributes to 
the m = 1 structure (see also Table HI and we assume the 
Kepler motion with radius « 10A/bh- The amplitude of 



the enhanced peaks could be larger than the noise level 
of the advanced ground-based detectors. For the SMBH 
formation scenario with (Mbhi^*) — (10^AfQ,10 Gpc), 
the predicted peak amplitude is h^s ^ 10~"'^^-10~^^ with 
the frequency ~ 10"^ Hz. Such GWs would be detected 
with a high signal-to-noise ratio > 10 by LISA. 

Summary and discussion. — We have explored the PPI 
in BH-torus systems in the framework of numerical rel- 
ativity. We have found that the to = 1 PPI occurs for 
a wide range of self-gravitating tori orbiting BHs. The 
resulting m — 1 structure in the torus survives for a long 
time after the PPI growth saturation, and large ampli- 
tude, quasiperiodic GWs are emitted. For reasonable hy- 
pothetical masses and distances, we have found that the 
quasiperiodic GWs emitted by this mechanism could be 
a promising source for ground-based and spacecraft de- 
tectors (for stellar-mass BHs and SMBHs, respectively). 

The BH-torus system composed of a stellar-mass BH 
is a promising candidate for the central engine of GRBs, 
and our results suggest that the so-called collapsar hy- 
pothesis 6] may be verified via observation of GWs. In 
GRBs, jets are likely to be ejected along the rotation 
axis of a spinning BH and torus. If the BH-torus sys- 
tem formed in the collapsar is axially symmetric, GWs 
of Z = 2 and to = mode are primarily emitted in the di- 
rection perpendicular to the rotation axis (e.g. [lij). By 
contrast, the quasiperiodic GWs studied in this letter are 
emitted along the rotation axis because the I ^ m = 2 
mode is dominant. Therefore, it could be possible to 
observe GRBs and GWs simultaneously, and to explore 
the collapsar scenario via GW observation. On the other 
hand, in the formation of a SMBH following SMS col- 
lapse, the torus-to-BH mass ratio is predicted to be 0.05- 
0.1 2\. Our results suggest that the PPI could grow in 
such a system. This implies that this scenario may be 
confirmed by the GW detector LISA. 
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